Numerical simulations on shear behaviour of rock joint network under constant normal stiffness conditions

In this study, the numerical direct shear tests were conducted to investigate the shear mechanical properties of joint networks under constant normal stiffness (CNS) boundary conditions. The influence of random joint number on shear stress (τ), dilation (normal displacement, δv) and normal stress (σn) of rock mass were studied quantitatively with fixed main slip surface. At the same time, the internal stress evolution process and failure process were analyzed. The results reveal that the number of random joints (γ) has little effect on the shear and normal stresses. The normal displacement of the sample generally decreases as the number of random joints increases. In addition, the normal displacement of the specimen is absorbed by the random joints when the number of random joints in the specimen increases to a certain level: when γ is greater than 6 and the shear displacement (μ) reaching 10 mm, the specimen exhibits shear contraction. Therefore, the internal random joints mainly control the failure mode and dilatancy performance of the specimen, while the main joint of the rock controls the shear stress of the specimen.


Introduction
Shear failure is a typical failure form of the surrounding rock in high-stress underground engineering [1]. The shear cracks developed in the surrounding rock of high-stress underground engineering cut the complete surrounding rock into blocks [2]. Under the action of the redistributed stress of the surrounding rock, these blocks undergo shear slippage along the fracture surface, resulting in disasters, such as large deformation of the roadway and damage to the supporting structure [3]. The shear performance of a rock mass is strongly affected by the presence of joints [4]. Recognizing the significance of its evaluation for underground engineering, many researchers have focused on the shear characteristics of the joint surfaces since the 1960s. For deep underground rock engineering, when a shear-slip failure occurs, the normal stress acting on the rock mass would continue to increase due to the restriction of the surrounding rock in the underground or rock bolts and cannot freely shear dilatation [5][6][7]. Therefore, constant normal load (CNL) boundary conditions can no longer meet the site requirements [8], and constant normal stiffness (CNS) boundary conditions that are more in line with the actual situation should be adopted [9]. Therefore, it is of great significance to investigate the shear behaviour of joints under CNS boundary conditions by designing shear tests that conform to the stress state of surrounding rock in deep underground rock engineering [10].
To better investigate the direct shear characteristics of jointed rock mass, some scholars have carried out experimental researches. Indraratna et al. [11,12] carried out research on the shear behavior of infilled joints under the boundary conditions of CNS based on the world's first independently developed CNS shear equipment. The research results showed that even if a thin layer of filler was added, the shear strength of the joint would be significantly reduced. Since the equipment of Indraratna et al. uses spring to simulate the normal stiffness, their equipment cannot realize the free switching of the normal stiffness. Subsequently, Jiang et al. [13,14] developed the first servo-controlled CNC shearing machine, which realized the free switching of boundary conditions. Jiang et al. found that the normal stiffness has a significant influence on the mechanical behavior of joints during the shear process based on the shear tests conducted using the servo-controlled CNC shearing machine. Xia et al. [15] developed a shear-flow coupling experiment system for rock joints under CNS boundary conditions, and studied the hydraulic coupling characteristics in the shear process. In addition, the shear behavior of joints under CNS boundary conditions has been further studied and developed with the support of the above-mentioned equipment [16][17][18][19].
In summary, previous studies have reflected the shear mechanical properties of rock joints under CNS boundary conditions, including shear strength, dilatancy characteristics and joint surface damage. However, most of the studies on rock joint shear under CNS boundary conditions [20] are limited to the study of single joints due to the test conditions and difficulties in specimen preparation. But joints do not exist in isolation in rock masses, and they are often widely distributed within the rock masses in the form of networks [21]. Since indoor experiments have great difficulties in preparing joint network rock masses, numerical simulation studies are needed to fill the gaps in this research field. The numerical simulation method (the Finite Element Method, FEM or the Distinct Element Method, DEM) [22], laboratory test and theoretical analysis are known as the three main analysis methods. Particularly, the repeatability, simplicity, convenience and diversity of numerical simulation methods play an important role in the field of studying deep underground engineering [23]. For example, Xu [24] and Saadat M [25] used the particle flow code (PFC) to study sawtooth filling joints with different normal stresses, different undulation angles and different filler thicknesses. And the shear mechanical properties and crack propagation phenomena of the sawtooth filling joints were reproduced. Asadi M S [26,27], Bahaaddini [28,29] and Gutiérrez-Ch J.G [30] investigated the mechanical properties of shear damage and the evolution of micro-cracks during the shear process of serrated joints with different geometrical characteristics at different initial normal stresses (σ n0 ), different shear directions by using PFC. However, there are few numerical simulation studies on the shear properties of joint network under CNS conditions. Therefore, we adopted an innovative cyclic assignment method of the smooth joint model (SJM) [31] in this study to overcome the failure problem of the joint model during the shear test. The numerical simulation model of the jointed rock mass under CNS boundary conditions was firstly established to investigate the shear behavior of the joint network under the CNS boundary condition. Then, the influence of random joint number on shear stress, normal stress and normal displacement of jointed rock mass were analyzed. Grants (LQ21E040003), National Natural Science Foundation of China under Grants (52204146) and China Postdoctoral Science Foundation under Grants (2022M712301). " By the way, Yuan Gao is recipient of the first three fundings (under Grants: 22KIB560010, JC12022098, PCMGH202201) and he takes a role in data analysis and preparation of the manuscript; Guansheng Han is recipient of the fourth funding (under Grants: LQ21E040003) and he takes a role in study design, data collection and analysis, decision to publish and preparation of the manuscript; Yu Zhou is recipient of the last two fundings (under Grants: 52204146, 2022M712301)and he takes a role in data collection and analysis of the manuscript. Contributions of other authors can be found at the end of the manuscript.

Establishment of rock single joint model
In previous studies, the numerical simulation research of joint shear mainly adopted the following two methods: 1) Bond deletion method. This method represents the joint surface by removing the particle bonding within a certain distance on both sides of the joint trajectory. Although this method is simple, the resulting joint surface will cause joint roughness at the microscopic scale, resulting in interlocking particles on the joint surface will make the simulation results far from the actual results [32]. 2) The smooth joint model method. It is widely used in PFC to simulate the mechanical behavior of rock weak surface and the direct shear test. However, the direct application of the SJM cannot completely and correctly simulate the shear behavior of real joints. For example, Bahaaddini et al. [28,29] found that even a plane joint with a SJM can only correctly simulate the real shear strength when the shear displacement is smaller than the minimum particle radius. When the shear displacement is larger than the minimum particle radius, the simulated strength is greatly biased high or low. This is because the particles on both sides of the joint surface cannot be correctly identified during the shear process when the SJM is directly applied, and the resulting stress concentration in local particle contact will cause distortion of the shear stress. Based on the above deficiencies, Bahaaddini et al. adopted a new shear box generation method and an SJM assignment method, and achieved good application results [33][34][35][36][37]. In this work, the numerical simulation of joint shear under CNS boundary conditions is also studied based on this method.
Taking the example of a single joint (the contour curves of standard joint roughness coefficient (JRC) were proposed by Barton and Choubey [38], JRC = 10-12), the establishment steps of a single joint model in deep rock are as follows: Initial shear box assembly (see Fig 1A): Use the wall command in PFC to generate the upper and lower parts of the frictionless wall separately. At the same time, the geometry command is used to import complex joint contours. Since only one side of a wall is activated and effective, the generation of a complex joint profile requires the creation of two walls.
Generation of particles (see Fig 1B): Uniformly distributed particle aggregates satisfying a certain pore density are formed in the upper and lower walls respectively. The minimum radius (R min ) is 0.5 mm and the maximum radius (R max ) is 0.75 mm. Such particle size not only satisfies the calculation requirements of computer capacity but also has relatively accurate and stable calculation results. This step requires sufficient time to make the generated particle assembly reach a static equilibrium state.
Apply the Parallel Bond Model (see Fig 1C): In this step, the floating particles (particles with less than three contacts) in the particle assembly generated in step 'Generation of particles' are deleted first to improve computational efficiency. Then the parallel bond model (PBM) is used to simulate the mechanical behavior of the intact rock part. And the linear contact model (LM) was assigned between particles and walls.
Assignment Smooth Joint Model (see Fig 1D): The application of the SJM is the core of this modeling approach. Firstly, the joint surface between the up and down particle groups is identified. Then the smooth joint model is continuously assigned to the joint surface during the application of normal stress and shear, and the smooth joint direction at each contact is consistent with the overall joint. This approach can ensure that there must be a smooth joint model where the joint surfaces are in contact during the shear process. Otherwise, the effectiveness of the joint model would cause distortion of the simulated shear strength as the μ increases.
The numerical model and boundary conditions are shown in Fig 1A-1E. Corresponding calibrated micro-parameters are listed in Table 1.

Calibration of parameters
In the physical test, the super rock brand plaster produced in Japan is used as the rock-like materials. The basic mechanical properties of the plaster are as follows: the optimal water mixing ratio is 1:0.2; The hardening time is 10 minutes; the expansion rate is 0.09% in 2 hours; the expansion rate is 0.13% in 24 hours; the uniaxial compressive strength after 1 hour is 50 MPa. Therefore, the samples cast with plaster can better simulate the real rock situation.
First, PBM and LM were calibrated and validated by uniaxial compression tests (see Fig 2A  and 2B), the results show that the uniaxial compressive strength, peak strain, Young's modulus and Poisson's ratio of the complete specimen in the numerical simulation test are basically consistent with the physical model test. Then SJM was calibrated by CNS shear test (see Fig 2C  and 2D). When JRC = 10~12 and σ n0 is 1, 3 and 5 MPa respectively, the shear stress-shear displacement curve and the final failure mode of the specimen in the numerical simulation test are basically consistent with the physical test. Therefore, we believe that the micro-parameters adopted in the numerical simulation experiments are reasonable.

Establishment and numerical calculation scheme of rock random joint network model
The size of the rock random joint network model was designed as 200mm×100mm. In order to simplify the calculation model, the rock random joint network model mentioned in this Table 1. Micromechanical parameters of PFC numerical model.

Micro-parameters Values
Particle and linear contact Radius of minimum particle, R min (mm) 0.5 Particle-size ratio, R min / R max 1.5 Particle density, ρ (g/cm 3 ) 1546 Coefficient of particle friction, υ 0 Smooth-joint contact friction coefficient, υ sj 0.45 paper consists of a main joint and random joints. Standard JRC Curve with JRC = 10~12 for the main joint. The random joints are through joints, which are simulated by deleting particles at specific positions. The length is set to 50 mm, and their spatial positions are randomly generated and conform to the Gaussian function distribution. In addition, the number of random joints is the research variable of this paper. The number of random joints is expressed by γ, and a total of 12 working conditions (γ increases from 1 to 12) are studied. A partial rock random joint network models are shown in Fig 3. The specific test scheme is as follows: the initial normal stress is set to 1MPa, the normal stiffness (k n ) is set to 3 GPa/m, the maximum shear displacement (10 mm) is 5% of the specimen length (200 mm) and the shear rate is set to 0.1 m/s. (The loading speed of 0.1'm/s' can fully simulate the quasi-static loading conditions; when the loading speed is 0.1'm/s', it takes about 3750 steps to produce 1mm shear displacement).

Influence of random joint number on shear behaviours
In terms of shear stress: under CNS boundary conditions, the variations of shear stress of rock random joint network specimens under different random joint numbers are shown in Fig 4A.  It can be seen from Fig 4A that the number of random joints has little effect on the shear stress of the specimen. The cures of shear stress versus (vs.) shear displacement are mostly coincides while the number of random joints increases from 1 to 12. It can be seen that the main joint controls the shear stress of the sample to a large extent.
The curve of shear stress vs. shear displacement can be roughly divided into rapid rise stage, serrated decline stage and residual stage. Taking γ = 2 as an example (Fig 4B), when the μ increases from 0 to 1.16 mm, the τ of the sample reaches the first peak value of 1.59 MPa (Stage I). When the μ is between 1.16 mm and 5.6 mm, the shear stress of the specimen appears multiple peaks successively (all of which are lower than the first peak in value), and the peak intensity gradually decreases (Stage II). When the μ increases from 5.6 mm to 10 mm, the shear stress of the specimen fluctuates within a range of fewer than 0.53 MPa (Stage III).
On the one hand, the shear stress of the specimen decreases due to the shearing of the convex joint surface during the shear test. On the other hand, the shear stress of the specimen  increases due to the constraint of normal stiffness. The two factors compete with each other, so the phenomenon of the serrated decrease in shear stress appears.
In terms of normal stress: the evolutions of σ n of rock joint network specimens under different number of random joints under CNS boundary conditions are shown in Fig 5A. It can be seen from Fig 5A that the number of random joints has little effect on the normal stress of the specimen. The trends of the curve of normal stress vs. shear displacement are basically consistent when the number of random joints increases from 1 to 12.
With the increase of μ, the normal stress of the specimen can be divided into rapid rise stage, stable stage and residual stage. Also taking γ = 2 as an example (Fig 5B): Stage A: When the μ increases from 0 to 3.3 mm, the normal stress of the specimen increases rapidly from 1 Ma (σ n0 ) to the first peak value of 2.7 MPa. Stage B: When the μ increases from 3.3 to 5.6 mm, the normal stress of the specimen fluctuates in the range of 2.2 MPa to 2.7 MPa. Stage C: When the μ increases from 5.6 to 10 mm, the normal stress of the specimen fluctuates in the range of 1.5 MPa to 2.1 MPa. In the CNS shear test, the normal stress of the specimen is generally higher than the initial normal stresses due to the shear dilatancy of joints (Eq 1) [13].
Where σ n is the real-time normal stress of the specimen, σ n0 is the initial normal stress set by the test, k n is the normal stiffness set by the test, and δ v is the normal displacement increment.
In terms of normal displacement: the variations of normal displacement of rock joint network specimens under different numbers of random joints under CNS boundary conditions are shown in Fig 6A. With the increase of μ, the normal displacement of the specimen can be divided into rapid increase and gradual decrease stages. As the number of random joints increases, the normal displacement of the specimen decreases gradually. It is noteworthy that the specimen exhibits shear contraction when the μ is equal to 10 mm when the number of random joints in the specimen is greater than 6. Taking γ = 1 and γ = 12 as examples (Fig 6B): during the whole shear test, the normal displacement of the specimen with γ = 1 is higher than the specimen with γ = 12; When the μ is equal to 10 mm, the final state of the sample with γ = 1 shows shear dilatancy, and the dilatancy amplitude is +1.06×10-4 m. The final state of the sample with γ = 12 shows shear contraction, and the shear contraction amplitude is -1.70×10-5 m. The above phenomenon shows that after the number of random joints in the specimen reaches a certain extent, the dilatancy generated during the shear process of the specimen will be absorbed by the random joints.

Micromechanical behavior: Evolution of inner stress and the failure modes
The inner stress evolution and failure mode of the specimen during the shear process are shown in Figs 7 and 8. Figs 7 and 8 are the inner stress evolution process and failure process of the rock joint network specimens corresponding to γ = 1 and γ = 12. Among them, (a)~(f) are the stress distribution of the specimen at the μ is 0.2, 2, 4, 6, 8 and 10 mm, respectively. The red area indicates the compression area, and the blue area indicates the crack. And (g) is the relationship between shear stress, tensile crack and shear displacement.
The inner stress distribution of the specimen corresponding to γ = 1 is shown in Fig 7. When the μ is small (μ = 0.2 mm), the interior of the specimen is compressed. There is no obvious stress concentration area and obvious microcrack. When the μ increases to 2 mm,  microcracks appear at the protruding position of the joint, and the stress concentration area appears nearby. With the increase of μ, the internal compression area of the specimen decreases gradually and the number of microcracks increases significantly. The specimen is gradually destroyed. When the μ reaches the maximum value (10 mm), the failure area of the specimen is concentrated near the main joint. When γ = 1, the existence of random joints has little effect on the failure mode of the specimen. In addition, the microcracks inside the specimen are mainly tensile cracks.
The inner stress distribution of the specimen corresponding to γ = 12 is shown in Fig 8. When the μ is equal to 2 mm, microcracks are generated inside the specimen at the protrusion position of the main joint and the tip of the random joint. With the increase of μ, the number of microcrack areas is increasing, and the number of microcracks is also increasing. It is worth noting that the locations of cracks are near the main joint or the tip of the random joint. When the μ is equal to 3 mm, the τ of the specimen begins to decrease obviously, and the number of tensile cracks inside the specimen increases rapidly.
Meanwhile, the shear resistance of the specimen is affected by the rapid development and penetration of the microcracks. When the μ is equal to 6 mm, the specimen almost loses the shear resistance. When the μ continues to increase, the number of tensile cracks in the specimen increases slowly. Some tensile cracks will be generated at the tips of random joints, and these tensile cracks will continue to develop and penetrate during the shear test. The existence of random joints affects the failure mode of the specimen to a certain extent.

Conclusions
The shear effect of rock random joint network specimens under CNS boundary conditions was studied by numerical simulation. The reproduce of CNS boundary condition, the rough joints and the cyclic assignment of joints have been well realized through the self-developed PFC2D code. The mainly obtained conclusions are as follows: 1. The number of random joints has little effect on the shear stress and normal stress of the sample. The variation curve of shear stress with shear displacement can be roughly divided into rapid rise stage, serrated decline stage and residual stage. With the increase of μ, the normal stress of the specimen can be divided into rapid rise stage, stable stage and residual stage. As the number of random joints increases, the normal displacement of the specimen generally decreases.
2. The shear dilatancy generated during the shear process of the specimen will be absorbed by the random joints after the number of random joints in the specimen reaches a certain extent. For example, when the δ h is 10 mm, the final state of the sample with γ = 1 shows shear dilatancy, and the shear dilatancy amplitude is +1.06×10-4 m. And the final state of the sample with γ = 12 shows shear contraction, and the shear contraction amplitude is -1.70×10-5 m.
3. In the shear test of rock random joint network, the random joints inside the rock mainly affect the failure mode and the shear dilatancy performance of the specimen, while the main joint of the rock controls the shear stress of the specimen.